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We present a computer simulation study of glassy and crystalline states using the standard 
Lennard-Jones interaction potential that is truncated at a finite cut-off distance, as is typical of 
many computer simulations. We demonstrate that the discontinuity at the cut-off distance in the 
first derivative of the potential (corresponding to the interparticle force) leads to the appearance 
of cut-off nonlinearities. These cut-off nonlinearities persist into the very-low-temperature regime 
thereby affecting low-temperature thermal vibrations, which leads to a breakdown of the harmonic 
approximation for many eigen modes, particularly for low-frequency vibrational modes. Further¬ 
more, while expansion nonlinearities which are due to higher order terms in the Taylor expansion 
of the interaction potential that are usually ignored at low temperatures and show up as the tem¬ 
perature increases, cut-off nonlinearities can become most significant at the lowest temperatures. 
Anharmonic effects readily show up in the elastic moduli which not only depend on the eigen frequen¬ 
cies, but are crucially sensitive to the eigen vectors of the normal modes. Whereas, those observables 
that rely mainly on static structural information or just the eigen frequencies, such as the vibra¬ 
tional density of states, total potential energy, and specific heat, show negligible dependence on the 
presence of the cut-off. Similar aspects of nonlinear behavior have recently been reported in model 
granular materials, where the constituent particles interact through, finite-range, purely-repulsive, 
potentials. These nonlinearities have been ascribed to the nature of the sudden cut-off at contact in 
the force-law, thus we demonstrate that cut-off nonlinearities emerge as a general feature of ordered 
and disordered solid state systems interacting through truncated potentials. 

PACS numbers: 02.70.-c, 63.50.-x, 63.20.-e, 62.25.-g 


I. INTRODUCTION 


Let us consider a Hamiltonian system, composed of 
N point particles, in a 3-dimensional solid (crystalline 
or glassy) state. We assume that the potential energy, 
$ = $({u}), of the system is a function only of the 
positions of the particles, {^ 1 , 7 ^ 2 , ...,r)v}. When the am¬ 
plitude of vibrational motions of the particles is small 
enough, the harmonic approximation provides a suitable 
description of the system [l| , where the equations of mo¬ 
tion are linearized 




(7 = l,2,...,iV). (1) 


Here, rm is the mass of particle i, Uiit) = ri{t) — is the 
displacement of particle i from the potential minimum 
state “0” where <f> = $° = $({ 7 =) = }), and ()° denotes 

the value at the minimum. The solution of Eq. 0 can 
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be written as the superposition of 3(iV — 1) eigen modes 
(3 zero-frequency translational modes are removed) 

3N-3 

= '^ A^exp{juj^t)^ {i = l,2,...,N), (2) 

fe=i 

where Aq is the amplitude of mode k, and j is the imagi¬ 
nary unit, and ^ = [e^, e^,..., are respectively the 
eigen frequency and the eigen vector of mode k, with 
the polarization vector of particle i, which are obtained 
by solving the eigenvalue problem of the 3N x 3N dynam¬ 
ical matrix, or Hessian, /dridvjf' [i[. For thermal 

systems, e.g. crystals and glasses, in the low tempera¬ 
ture regime such that thermal motions of the particles 
are constrained within their local energy minima, the 
harmonic approximation successfully describes both vi¬ 
brational properties of the systems, and other thermody¬ 
namic observables, e.g. the specific heat and elastic con¬ 
stants. The success of this picture relies on the notion 
that particle vibrations are analogous to a ball-spring 
model whereby interacting particles remain in ‘contact’ 
in the sense that particles always feel the interactions of 
their neighbors. 

In contrast, recent studies of purely-repulsive, short- 
range interaction potentials as model systems of partic- 
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ulate matter have revealed a strong nonlinear or anhar- 
monic feature in their vibrational modes This non¬ 

linear behavior has been ascribed to the nature of the 
sudden cut-off at contact in the force-law between non- 
cohesive granular particles. Indeed, given the one-sided 
nature of such potentials and the fact that their second 
derivatives are non-analytic at the contact distance Q, 
once a contact altering event (opening or closing) occurs, 
the harmonic approximation, Eq. o, cannot c^ture this 
energy variation, and nonlinearities emerge (2l-l7|. We 
denote these aspects of the interaction law as “cut-off 
nonlinearities” and emphasize that these are quite dis¬ 
tinct from the usual “expansion nonlinearities” which 
arise from higher order terms in the Taylor expansion 
of the potential energy, $ = ‘f’dr)}), around the energy 
minimum 

In the present contribution, we show that cut-off non- 
linearities are not specific to granular matter, but appear 
as a general feature in any solid state system in which 
the interaction is cut off at a finite distance, as is com¬ 
mon practice in computer simulations. While enforcing 
the truncation of quasi-long-ranged interaction poten¬ 
tials, e.g., attractive Lennard-Jones potentials, is largely 
a numerical issue to improve computational efficiency, 
from a mathematical perspective the appearance of cut¬ 
off nonlinearities that results from this truncation share 
the same mechanism with those in particulate, jammed 
systems. Therefore, it seems pertinent to ask whether 
such cut-off nonlinearities actually affect structural and 
mechanical properties of solid state systems. We will 
show that, although simple structural properties are not 
influenced by these nonlinearities, they become appar¬ 
ent in normal modes analyses. In particular, such cut-off 
nonlinearities can strongly affect mechanical properties 
that depend explicitly on the vibrational modes, specif¬ 
ically, the eigen vectors and thence the particle polar¬ 
ization vectors. Oddly, this effect can be enhanced by 
lowering the temperature. 

We study the Lennard-Jones (LJ) pair potential 
’pLji'f) (our results for the LJ potential can be extended 
to more practical potentials, and should be quite general) 


^Ljir) = 4e 

= - 


Y-V- 

(-Y 

W / 

\ r / 


24e 

r 


(7)"-(7)' 


( 3 ) 


where ' denotes the derivative with respect to distance 
r. The values of e and a characterize the energy and 
length scales of the interaction, respectively. The stan¬ 
dard protocol in numerical simulations is to truncate the 
potential at some cut-off distance r = m 

<i)TLj{r) = (j)Lj{r)H{rc-r), 

<t>TLAr) = (l>Lj{r)H{r^ - r) - (t>Lj{r)S{r - r^), 

where H{x) is the Heaviside function {H{x) = 1 for 
a; > 0, Hix) = 0 otherwise), and 5{x) is the impul¬ 
sive function. This truncation leads to a discontinuity 


at r = Tc in the potential (jiTLjiA^ impulsive 

term, ~ S{r — rc), in its first derivative i.e. the 

interparticle force. To prevent the discontinuity in the 
potential and impulsive force, it is standard practice to 
use the “shifted potential” (/>sp(r') Eini; 

<l>sp{r) = [ 4 >Lj{r) - (l>Lj{rc)\ H{rc - r), 

AspA) = ALj{r)H{r,-T). ^ > 

The function (j)sp{r) is now continuous at r = Tc and 
has no impulsive term in its first derivative, but still has 
discontinuities at r = Tc in its derivatives, including the 
force term. To smooth the potential further, e.g. make 
the first derivate continuous, we can employ the “shifted- 
force” potential ^!>sf(?’) Elll; 

(I^SFir) = AljA) - (j)Lj{rc) -{r - rc)(t>'Lj{rc)\H{rc - r), 

Asp A) = WljA) - ALj{rc)\H{r^ - r). 

( 6 ) 

We can also smooth arbitrarily high nth-order (n > 2) 
derivatives by adding terms, ^ (r — Tc)", in the same 
way as in (pspA)- Another option is to interpolate the 
potential around the cut-off distance r = rc by using 
polynomials (see e.g. Ref. M)- 

If the second derivative of the truncated potential 
is non-analytic at r = as it is the case for all of 
4>tlj{i’)Asp{t)Asp(t’)^ if some pairs of particles 
pass through r = Vc, the dynamical matrix and Eq. (HD 
cannot capture the energy change due to this event, and 
as a result, cut-off nonlinearities emerge. In the following, 
we will show that the cut-off discontinuity in the inter¬ 
particle force (the first derivative of the potential) en¬ 
hances the cut-off nonlinearities, and has non-negligible 
effects on the low-temperature thermal vibrations, even 
though the discontinuity is small. Specifically, when we 
use the shifted-potential (()sp(r), where the interparticle 
force is discontinuous at r = rc, the harmonic approxi¬ 
mation breaks down even at very low temperatures for 
many vibrational eigen modes, particularly those that lie 
at lower frequencies. However, when the shifted-force 
potential (j^spA) i® employed to make the interparticle 
force continuous, then the cut-off nonlinearities are sig¬ 
nificantly suppressed, and the harmonic approximation 
becomes applicable again. We also study the effects of 
cut-off nonlinearities on several structural and thermo¬ 
dynamic observables. In particular, we will show that 
the cut-off nonlinearities have a negligible effect on those 
properties that do not depend on the details of the eigen 
vectors of the normal mode decomposition. These in¬ 
clude the radial distribution function (RDF), vibrational 
density of states (vDOS), total potential energy, and spe¬ 
cific heat. Whereas, mechanical properties such as the 
elastic constants, which depend crucially on particle po¬ 
larization vectors, are much more strongly affected. 

The paper is organized as follows. In Sec. HIl we de¬ 
scribe in detail the numerical systems studied. SectionHIIl 
contains a comprehensive presentation of our results. We 
first discuss the results of the static structure and vibra¬ 
tional states in the harmonic limit T = 0 in Sec. lHI Al We 
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Potential Force at r = rc System /N 

1- (pspir) with Tc = 2.5 Discontinuous Glassl —6.85281 

Crystall —7.36841 

2. 0SF(r) with Cc = 2.5 Continuous Glass2 —6.07778 
_Crystal2 -6.60853 

3. 4>sp{''') with Cc = 3.0 Discontinuous Glass3 —7.29178 

Crystal3 -7.81922 


TABLE I. The potentials and the systems investigated in this 
work. We consider three types of truncated LJ potentials, 
1, 2, and 3, and two types of configurations, glass and face 
centered cubic (FCC) crystal, therefore a total of 6 = 3 x 2 
systems, which are listed in the table. We also present the 
value of the minimum potential energy per particle, 4>°/A^, 
for each system. The detailed descriptions of the potentials 
and the systems are found in the main text. See also Fig. [T] 


next study the effect of cut-off nonlinearities on thermal 
vibrations at finite temperatures T > 0 in Sec. IIIIBl and 
their effects on several physical quantities in Sec. IIII Cl 
We summarize our results in Sec IlYl and give our con¬ 
clusion in Sec. El 


II. SYSTEMS DESCRIPTION 
A. Lennard-Jones potentials 

We consider three types of truncated LJ potentials, as 
listed in Table |T] and plotted in Fig. [T] Previous numer¬ 
ical works usually made the potential (t>(r) continuous 
(shifted) at the cut-off distance r = Pc [TO - fo , in 

order to avoid the impulsive force, ^ 6(r — Pc). Then, 
each work optionally took care of discontinuities in the 
derivatives of the potential. Even if the potential (/)(p) 
was not made continuous (not shifted) at r = r^, the 
impulsive force was usually ignored, although the inter¬ 
particle force and ^^(r) become inconsistent in this case. 
Therefore, the shifted potential 0sp(r) (see Eq. (O) is 
considered to be the most practical to implement and is 
regarded as the standard potential to use in simulation 
studies Our focus in this study is then 

the discontinuity in the first derivative of the interac¬ 
tion potential (interparticle force), which in turn has a 
non-negligible effect on the low-temperature vibrational 
motions of particles. 

Potential 1: As the first potential “1” we use the 
shifted potential, ^sp(r), with Pc = 2.5. This potential 
has a discontinuity in its first derivative at r = = 

2.5. The value of the cut-off distance, Pc = 2.5, has 
been employed in many previous simulation works [ 13 - 
fl3| . Note that the value of original potential (/>lj(p) at 
Pc = 2.5 is 1.6% of the energy scale e. Our previous 
simulations of LJ glasses also used this potential 1 (with 
Pc = 2.5) to study acoustic excitations da and spatial 
distributions of local elastic moduli [l3| • 



1 1.5 2 2.5 3 3.5 

r 

FIG. 1. (a) Lennard-Jones potential (/)(r) and (b) its first 

derivative (/>'(p)- In the present study, we use three types of 
truncated potentials; 1: cf>sp{p) with Pc = 2.5, 2: cf)sF{r) with 
Pc ~ 2.5, and 3: 4>sp{p) with Pc = 3.0, as listed in Table U 
(see also Eqs. ([U for 4>sp{p) and (|6]) for (f>sF{p))- The red 
solid, green dashed, and blue dotted lines show the potentials 
1, 2, and 3, respectively. The insets show close-ups around 
the cut-off distance r = r^. 


Potential 2: In order to consider the continuous in¬ 
terparticle force, we use, as the second potential “2”, the 
shifted-force potential, ^ 5 p(r), with the same cut-off dis¬ 
tance Pc = 2.5 (see Eq. (HD- 
Potential 3: For our third potential “3”, we consider 
a longer cut-off distance Pc = 3.0 in the shifted potential 
4'spij')- The interparticle force has a discontinuity at 
p = Pc, but it is reduced by the longer Pc = 3.0, compared 
to that of potential 1 with Pc = 2.5, as we can visually 
recognize in the inset of Fig. [TJb). 

By comparing the results from those three potentials, 
we study the impact of the discontinuity in the inter¬ 
particle force at r = Pc, on cut-off nonlinearities and 
finite-temperature thermal vibrations. 


B. Systems preparation 

For each LJ potential described above, we prepared 
two types of configurations, an amorphous glass and 
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a face centered cubic {FCC) crystal, under periodic 
boundary conditions in the x, y, and z directions. The 
system is mono-disperse, composed of = 4,000 identi¬ 
cal particles with mass m, diameter a, and interparticle 
potential energy e. Throughout this paper, we present 
the values of quantities in units of a (length), e (energy), 
and m (mass). The temperature T and frequency w are 
expressed in units of e/ks {ks is Boltzmann’s constant) 
and T~^ = (ma^ j^ respectively. We fixed the num¬ 
ber density at p = NjV = 1.015, where the system length 
is L = = 15.8, and the lattice constant of the FCC 

crystal is a = L/10 = 1.58. The melting temperature Tm 
and the glass transition temperature Tg of this mono- 
disperse system have been reported as Tm — 1-0 and 
Tg ~ 0.4 O. 

We first considered potential 1, to prepare “glass 1” 
and “crystal 1” at a temperature T = 10“^, well below 
both Tg and T^. For the preparation of glass 1, we equili¬ 
brated the system at T = 2.0 in the normal liquid phase, 
and quenched it down to T = 10“^ with an extremely 
fast rate dT/dt = 4 x 10^, followed by an equilibration 
run at T = 10“^, as described in Ref. [l3|. After prepar¬ 
ing the glass and crystal with potential 1, we switched 
the potential from 1 to 2 and 3, and equilibrated the 
systems again at the same temperature T = 10“^. For 
the glassy cases, we relaxed the systems for sufhciently 
long times to eliminate any aging effects in these three 
cases. We note that all systems remain very close to the 
same energy minimum during the entire trajectory. After 
the systems were equilibrated at T = 10“^, we quenched 
them using the steepest descent method, from T = 10“^ 
to 0, i.e. into the nearest energy minimum 0 (inherent 
structure). 

The values of the minimum potential energy per parti¬ 
cle, <I)°/A^, presented in TablelH are different for the three 
potentials, with a difference of about 15% between the 
potentials 2 and 3, for both glass and crystal. Since dif¬ 
ferent cut-offs give different Hamiltonians, system prop¬ 
erties generally depend on the cut-off nature. In fact, 
as it has been reported in Refs. [IMil, thermodynamic 
properties, including the melting point Tm (and possibly 
the glass transition point Tg), strongly depend on the 
cut-off treatment. However, as we will see in the RDF 
g{r) of Fig. [5] in Sec. IHI A 11 the three truncated poten¬ 
tials produced practically identical configurations (glass 
or crystal). This allows us to focus entirely on the role 
of the cut-off nonlinearities in the same statie structures. 
Note that we switched to potentials 2 and 3 only after 
the system was prepared with potential 1; the potentials 
2 and 3 were not used from the initial stage where we 
generated liquid configuration at T = 2.0. The reason 
why we employed this particular preparation procedure 
was to minimize any possible differences in the configu¬ 
rations from which we started to compute observables. 

Finally, each system was heated from the energy min¬ 
imum state (T = 0). We study the low temperature 
regime, ranging from T = 10“^ to 10“^ for glasses, 
and from T = 10“^ to 10“^ for crystals. This tem- 




FIG. 2. The RDF g{r) in the energy minimum state (T = 0) 
for (a) glasses and (b) crystals. The values for three types of 
truncated potentials, 1 (red solid line), 2 (green dashed line), 
and 3 (blue dotted line), are presented. The inset to (b) is a 
close-up around r = 2.5. 

perature range is well below (one order of magnitude 
lower than) the glass transition Tg ~ 0.4 and the melting 
Tm — 1.0 temperatures. At each T, we carried out an 
equilibration followed by a production run, using NVT 
molecular-dynamics (MD) simulation. Here, we set the 
time step to 6t = 10“^ for T < 10“^, and 5t = 5 x 10“^ 
for T > 10“^. During the simulations, we observed no 
particle rearrangements, indicating that particles vibrate 
around the same energy minimum state at all studied 
temperatures, i.e. each system remains in its original 
inherent structure. All MD simulations were performed 
using LAMMPS [13,[HI. 


III. RESULTS 

A. Static structure and vibrational states in the 
harmonic limit, T = 0 

Before studying thermal vibrations at finite tempera¬ 
tures T > 0, we first look at the static structure and the 
vibrational states in the energy minimum state T = 0 
(the harmonic limit). We will see that the three trun- 





























cated LJ potentials produce practically identical static 
structures and vibrational states at T = 0. 
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1. Static structure 


Figure [2] presents the RDF, g{r), in the energy mini¬ 
mum state 0 (T = 0). As in Fig. [H^a), the three poten¬ 
tials show no differences in g{r) for the glasses. In the 
crystals, we observe slight differences between potential 
1 and the other potentials 2,3, as highlighted in the inset 
of Fig. mb). While crystals 2 and 3 have a single peak 
around each lattice distance, double peaks are apparent 
in crystal 1. In the FCC lattice structure with lattice 
constant a = 1.58, all the particles have their nearest 
neighbors located at a distance r = ( 1 / 10 / 2 ) x a ~ 2.5. 
In this situation, if the interparticle force has a discon¬ 
tinuity at r = Tc = 2.5, particles cannot be stabilized 
at the exact lattice positions, rather they are slightly 
displaced due to the discontinuous force, leading to the 
double peaks in g(r) for crystal 1. This effect also exists 
in glasses 1 and 3, i.e. some pairs of particles located 
at the cut-off distance r = Tc are not stabilized precisely 
at r = re- However, this feature is not noticeable in 
the glass due to the smoothing effect of the amorphous 
structure. 

Although the locations of particles in crystal 1 do not 
lie precisely at the lattice positions, they are approxi¬ 
mately located at these lattice points, with deviations 
of only 1%, hence they show the same FCC structure 
as crystals 2, 3. Therefore, we conclude that the three 
truncated potentials produce essentially identical static 
structures in the harmonic limit T = 0. 


2. Vibrational states 


As we explained in Eqs. © and harmonic vi¬ 
brational motions of particles are completely described 
in terms of the 3iV — 3 eigen frequencies u)^ and ei gen 
vectors e* = [^, ..., e^] (fc = 1,2,..., 3iV — 3) [l|. 

We diagonalized the Hessian {3N x 3N dynamical ma¬ 
trix), (9^d>/9r)9/)°, to obtain the eigen values and 
eigen vectors The eigen frequencies are given as 
and the eigen vectors are normalized such 
that e* • / = ■ ^i) = where Ski is the Kro- 

necker delta function. Note that if the system is stable 
in the energy minimum state 0, all the eigen values 
are positive, and the eigen frequencies are real num¬ 
bers. In this study, we used the ARPACK software [ 2 ^ 
to realize the diagonalization of the 3N x 3N matrix. 

Figure 13] shows the vDOS, d{uj), obtained from 3N — 3 
values of 


, 3Af-3 

= 3(jV-l) 




0 5 10 15 20 25 

UJ 

FIG. 3. The vDOS, d{uj) and d{uj), for (a) glasses and 
(b) crystals. The symbols represent the vDOS d(ca) obtained 
from the dynamical matrix (see Eq. ([7])), for potentials, 1 (red 
circle), 2 (green square), and 3 (blue triangle). The lines show 
the Fourier transform, d{u}), of the velocity auto-correlation 
function, d{t), at T = 10“® (see Eq. (I21II '). for potentials, 
1 (red solid line), 2 (green dashed line), and 3 (blue dotted 
line). The black solid line is the Debye prediction, which 
is calculated from the elastic moduli of potential 2. Note 
that the Debye prediction is almost unchanged for different 
potentials. 


In Fig. m we plot the participation ratio, PR^, for each 
mode fc, as a function of ui^; 


PR^ 


1 

/v 


■ N 

E'=*' 
. 2=1 


(et-et) 


-1 


( 8 ) 


which measures the extent of vibrational localization in 
mode k [22l - [25l| . The vDOS d(w) shows the typical shapes 
for glasses and crystals. Crystals: (i) at intermediate- 
high frequencies d(a;) has two branches, corresponding 
to transverse and longitudinal acoustic phonons, and (ii) 
at low frequency, Debye scalin g, d jui) ~ is observed 
(black solid line in Fig. |3Db)) (ll.l27j|. Glasses: (i) smooth 
variations and broadened distributions are observed, and 
(ii) the low frequency portion of d(uj) shows an enhance¬ 
ment of modes over the Debye prediction [l3, [13, ilslj . 
In addition, from Pi?^, we recognize localization of the 
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FIG. 4. Participation ratio PR^ versus eigen frequency 
for (a) glasses and (b) crystals. The values for the three po¬ 
tentials, 1 (red circle), 2 (green square), and 3 (blue triangle), 
are presented. 


low and high modes in the glasses, which are a 
characteristic feature of various amorphous [H, and 
jammed [2^ matters, whereas all the modes in the 
crystals are extended phonon modes. 

In Figs. [Handm there does not appear to be any no¬ 
ticeable differences among the three potentials. This is 
further validated by the results shown in Fig. |TT] (see 
lines) and related discussion. Therefore, the three trun¬ 
cated potentials produce practically identical vibrational 
states in the harmonic limit T = 0. We note, however, 
that small detailed differences are observed in c?(a;) and 
PR^] for instance, glass 2 with continuous interparticle 
force contains a number of lower frequency modes than 
those in glass 1 (and 3) with the discontinuous force (the 
lowest eigen frequency of = 0.79 in glass 2 is lower 
than that of uj^ = 1.13 in glass 1). These small differ¬ 
ences are picked up by the (non-afhne) elastic moduli 
which shows differences between the three potentials, as 
indicated in Fig. [15] (see lines). 

We close this section on the T = 0 properties with 
a brief discussion of the low frequency quasi-localized 
modes in glassy configurations, which will be useful in the 



FIG. 5. The spatial correlation function C^{rij) of eigen 
vector for several different eigen modes in glass 1. Note 

that the mode with = 1.13 and PR'^ = 0.11 is a low 
frequency quasi-localized mode, while the one with uj^ = 30.0 
and PR'^ — 0.08 is a high frequency localized mode [2^ - [^ . 


next section (Sec. IIIIB II) . We have studied the spatial 
correlation function C'^(ry) of the eigen vectors e*(r)) 
for each mode fc; 

C'=(r.,) = (et(r0.et(r)))(,,->, (9) 

where e* (r)) is a function of the position r), Vij = | r) — r} |, 
and denotes the average over all the pairs of parti¬ 
cles (ij). Figure [5] presents C^{rij) for several different 
eigen modes k in glass 1. It is observed that the lowest 
= 1.13 mode has the low value of PR^ = 0.11, but 
its C^{rij) shows a longer spatial correlation compared 
to other modes with higher values of PR^. Note that 
the negative correlation comes from the transverse na¬ 
ture of this mode. In the localized modes with low 
vibrational motions of particles are confined to a few 
small regions. Therefore one may expect that the spa¬ 
tial correlation between vibrational motions of particles 
is rather short within those confined regions. However, 
Fig. [5] confirms that the lowest = 1.13 mode actually 
shows an extended character with the long spatial corre¬ 
lation. Thus, the low mode is “quasi-localized” 

[ 2 ^ , which is distinct from “true” localization without any 
extended nature, observed in the high frequency modes. 
In Fig. IH the high frequency localized mode (w^ = 30.0 
and PR^ = 0.08) shows only a short spatial correlation, 
i.e. C^{rij) decays within ~ 1 (order of the parti¬ 
cle diameter). Similar observations have been found in 
athermal jammed solids [13, 


B. Cut-off nonlinearities in thermal vibrations at 
finite temperature T > 0 

We now turn our attention to thermal vibrations at 
finite temperatures T > 0. We focus on each eigen mode 
k and monitor vibrational motions of particles along the 
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mode k. Specifically, we measured the vibrational ampli¬ 
tude, A^[t), and the energy landscape, A$(yl^), for each 
mode fc, which will be presented in this section. When 
a pair of particles passes through the cut-off distance 
r = Tc, the cut-off nonlinearities affect the vibrational 
motions of particles. We will show that the discontinuity 
in the interparticle force at r = Tc enhances the cut-off 
nonlinearities, which causes a breakdown of the harmonic 
approximation for many eigen modes k. More precisely, 
the number of particles pairs that are involved in mode 
k and experience the force discontinuity determines the 
strength of the cut-off nonlinearities. 


1. Vibrational amplitude along eigen mode k 

We have extracted the vibrational amplitude, A^{t), 
along the eigen mode k from an MD trajectory in the 
(NVT)-ensemble, {ri (f), r 2 (f), ...,fAr(f)}, at each temper¬ 
ature T; 


N 

A'^{t)=J2Mt)-S^ (fc = l,2,...,3fV-3), (10) 

i=l 

where Ui{t) = fi{t) — r® is the displacement of particle i 
from its energy minimum position f*. We note that the 
MD trajectory is then described as the superposition of 
3A — 3 eigen modes k with the amplitude A^{t); 

3N-3 

u,{t)= Y. A'^it)^ (f = l,2,...,iV). (11) 

fc=i 


as follows. In the harmonic approximation, the poten¬ 
tial energy of the system. A®)!) = <l>(t) — measured 
from the minimum value is formulated in terms of the 
dynamical matrix, {d^^/dr'idrj)'^ 




N N 


i=l 

3N-3 


/ 92 $ 

V dridrj 


= E 


/c=l 

3N-3 


= E 


A^{t) 


[A>^{t)o 




EE 

i=i 3=1 
(,-i2 3N-3 


dndm 


■ e* 


— E 


fc=i 


k=l 


( 12 ) 

where Ep{t) = [A ^(t)Lc*^]^/2 is the potential energy of 
eigen mode k [291 - 131| . In the harmonic description, the 
total potential energy A$ is therefore the sum of Ep 
over all modes k. When the system is in thermodynamic 
equilibrium at a temperature T, energy equipartition 0 
implies that a potential energy T/2 is distributed equally 
among eigen modes, so that 


{El,{t)) =-T 


(AHtf) = 


T 




2 ’ 


(13) 


where, again, () denotes the NVT ensemble average 
(time average over the MD trajectory). Therefore, the 
harmonic value of the amplitude is given as 


A‘ 


harm 


■■=\ ) = 


Vt 




(14) 


In thermal systems at finite T > 0, the eigen modes al¬ 
ways exchange energy with each other, due to genuine 
expansion nonlinearities of the potential (and through 
the thermal bath as well) [I|, [29l - [^ . We remark that the 
expansion nonlinearities should be very small in our tem¬ 
perature range, but they still play a role in thermal equi¬ 
libration of the system. This fact is reflected in the obser¬ 
vation that all the amplitudes, A^{t) {k = 1,2,..., 3A—3), 
are a function of time t. As a consequence, the life-times 
of the modes, which can be measured by the mode energy 
correlation function [29l - l3lj| . are finite, although they be¬ 
come very large at low T. This mode mixing is a generic 
feature for T > 0 systems, and is active even at very 
low T, as in the present case. Particularly, in crystals, it 
is known as phonon-phonon interaction and plays a cen¬ 
tral role in thermal conduction (ll. [29l - [^ . In contrast, if 
Eqs. o and m hold, the modes can not interact. Thus 
the amplitudes A^{t) of modes are all constant, not de¬ 
pending on t, and the lifetimes of modes are infinite. Note 
that in Eq. © A>^{t) = A^exp{jujH), and \A^{t)\ = Ag 
is constant. Therefore, the strict definition of harmonic 
motions, Eqs. © and is never realized for systems 
at finite T > 0. 

Although strict harmonic vibrations are not realized in 
thermal systems, the law of equipartition of energy gives 
a harmonic condition for the amplitude A^{t) of mode k, 


where we explicitly denote the harmonic value by the 
subscript “harm”. Thus, we can identify harmonic vi¬ 
brations for mode fc, through the relation A^/= 1, 
where 



\j \ 




Vt/uj^ 


and deviations from the harmonic condition indicate the 
presence of anharmonicities. It is important to remark 
that at low temperatures, the expansion nonlinearities 
are very weak (even if they are still strong enough to 
induce mode interactions, as we mentioned above), and 
do not affect the harmonic approximation of the potential 
energy. Therefore, if for mode k, A^ /^ 1, this 
implies that the mode is deformed primarily by the cut¬ 
off nonlinearities. 

Figure |6] shows the ratio of the calculated amplitude 
and the harmonic expectation (solid line), A^/ 
versus the eigen frequency uj^ at the indicated values of T, 
for all investigated cases (glass and crystal). For poten¬ 
tial 1, many eigen modes k show larger amplitudes than 
the harmonic expectation, i.e. anharmonic vibrations. In 
glass 1, the low modes are anharmonic even at very 
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FIG. 6. Amplitude versus eigen frequency (a) glass 1, (b) glass 2, (c) glass 3, (d) crystal 1, (e) crystal 2, and (f) 
crystal 3. The value of is normalized by the harmonic expectation, = y/T/uj^, such that harmonic behavior of mode 

k corresponds to = 1 (Eq. (1151) 1. which is indicated by the horizontal solid line. The inset to (d) shows the large 

values of the amplitude in crystal 1. Different symbols represent three temperatures as indicated in the key of the middle panel. 


low T = 10“^, whereas at higher frequencies, uj^ > 5, the 
modes become harmonic. On the other hand, in crystal 
1 , many modes over the entire range of uj^ show large 
anharmonic effects. Remarkably, those anharmonic vi¬ 
brations are completely suppressed in potential 2, where 
all the modes are consistent with the harmonic expecta¬ 
tion value, in both glass 2 and crystal 2. We emphasize 
that cut-off nonlinearities should exist in both potentials 
1 and 2, since the second derivatives of the potentials are 
both non-analytic at r = Tc. What we have demonstrated 
here is that cut-off nonlinearities are most pronounced in 
the case of a discontinuity in the interparticle force at 
r = Tc, which indeed makes a visible impact on the ther¬ 
mal vibrations. As we will see below, the strength of 
cut-off nonlinearities in each mode k is then determined 
by the number of pairs of particles that are involved in 
the mode k and experience the force discontinuity. 

Comparing glass 1 and crystal 1: Crystal 1 shows 
much larger anharmonic amplitudes in many of the nor¬ 
mal modes. As we mentioned in Fig. [2jb) of the RDF 
g{r) (Sec. IIII A II) . all the particles in crystal 1 have 
their nearest neighbors located around the cut-off dis¬ 
tance r = Tc = 2.5 (which coincides with the exact 


lattice position), and in addition, the eigen modes are 
phonons, i.e. collective, extended modes. In this situa¬ 
tion, many modes have a large number of pairs of par¬ 
ticles which experience the discontinuous force causing 
large cut-off nonlinearities. On the other hand, in glass 
1 , fewer particles have such neighbors located precisely 
at r = Tc = 2.5 because of the amorphous structure. Fur¬ 
thermore, the eigen modes are less extended than regular 
phonons, exhibiting quasi-localized features, resulting in 
much smaller cut-off nonlinearities. 

Comparing glasses 1 and 3: We observe smaller 
anharmonic amplitudes in glass 3 than those in glass 1. 
In glass 3, although even more pairs of particles, located 
around = 3.0, are found than those at Tc = 2.5 in glass 
1 , the discontinuity in the interparticle force is smaller by 
the longer cut-off distance Vc = 3.0 (see Fig.[TJb)), which 
leads to the smaller cut-off nonlinearities. 

Comparing crystals 1 and 3: In crystal 3, all the 
modes show values around the harmonic expectation, i.e. 
the large anharmonic amplitudes observed in crystal 1 
completely disappear. The discontinuity at r = rc = 3.0 
in potential 3 is not located at the lattice points (see g{r) 
in Fig. [2Ib)), therefore no pairs of particles experience 
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FIG. 7. Mean square displacement of particles, (Ar^(t)) 
(Eq. (1161) 1. is plotted as a function of time t, for (a) glass 
1, (b) crystal 1, and (c) crystal 2. The value of (Ar^(t)) is 
normalized by the temperature T. Different lines represent 
three temperatures as indicated in the key. In the harmonic 
regime values of (Ar^(t)) /T for different T collapse onto a 
single curve. 


the discontinuity at low T s, and no cut-ofF nonlinearities 
appear in crystal 3. 

Low- modes in glasses 1 and 3: One noticeable 
observation in glasses 1 and 3 is that anharmonicities are 
found only in the low < 5 region, with lower modes 
showing larger anharmonic amplitudes, ^^/^harm > 
This result is explained by the fact that the lower 
modes involves more particles which experience the force 


discontinuity to cause cut-off nonlinearities. As we dis¬ 
cussed in Fig.[5]of (Sec.lIII A2I) . the low mode 

is actually quasi-localized [22 -1 ^ . i.e. although it is lo¬ 
calized, it exhibits an extended character with long spa¬ 
tial correlation. Therefore, the larger cut-off nonlineari¬ 
ties at the lower uj^ reflect the extended character of the 
low mode. 

Effect of temperature in crystal 1: Whereas 
at lower temperatures particle displacements are small 
enough that interacting pair separations remain close 
to the initial distance, at higher temperatures thermal 
vibrations increase in amplitude causing pairs of parti¬ 
cles, initially located near the cut-ofF distance r = rc, 
to pass through less often - particle separations reside 
inside or outside the cut-ofF distance - thereby reducing 
the effect of the cut-ofF nonlinearities. This is clearly 
observed in crystal 1, i.e. as T increases, anharmonic 
amplitudes are strongly suppressed toward the harmonic 
value, = 1. Here we note that expansion non- 

linearities are expected to be enhanced by larger vibra¬ 
tional amplitude at higher T, however we cannot resolve 
these effects at the values of T explored here. We also 
remark that cut-ofF nonlinearities have an effect that is 
opposite to the temperature of the expansion nonlinear¬ 
ities; increasing temperature leads to a suppression of 
cut-ofF nonlinearities, while expansion nonlinearities are 
expected to grow with temperature. 

Effect of temperature in glasses 1 and 3: The an¬ 
harmonicities in glasses 1 and 3 are quite insensitive to 
temperature. Given the fact that systems 1 and 3 share 
the same type of discontinuity in the force, one might ex¬ 
pect that the glassy systems should exhibit similar tem¬ 
perature effects as crystal 1, i.e. suppressed anharmonic¬ 
ities with increasing temperature. However, due to their 
amorphous structures, larger thermal vibrations also in¬ 
volve more pairs of particles that experience the force dis¬ 
continuity, and actually lead to an increase in the anhar¬ 
monic effects. These two effects - increasing temperature 
and a growing number of anharmonic vibrating particle 
pairs - appear to cancel, and as a result the anharmonic¬ 
ities in glasses 1 and 3 are less sensitive to temperature 
than might be expected. 

Finally, in Fig. [7] we plot the mean-squared- 
displacement (MSD) of the particles, 

(Ar^(t)) = -Y. W - • (16) 

In the harmonic limit, oc T. To highlight any 

deviations due to anharmonic behavior, we therefore plot 
the data normalized by T. We find excellent data collapse 
for both glass 1 (Fig. [7Da)) and crystal 2 (Fig. [TKc)). In 
glass 1, despite the fact that many low-frequency modes 
are anharmonic, these anharmonicities are, in a sense, 
hidden since all the 3N — 3 modes are mixed in (Af^(t)). 
On the other hand, we see clear anharmonic effect in 
(Af^(t)) for crystal 1 at the lowest T = 10“^ (Fig.EKb)); 
(Af^(t)) shows an increase of one order of magnitude 
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J^k^k 


FIG. 8. The potential energy landscape, A<1?(^^, j4o), 
along (a) low = 1.13 and (b) high u)^ = 8.79 modes in 
glass 1, and (c) low uj'‘ = 0.79 mode in glass 2. We stat¬ 
ically displace the particles from their minimum positions, 
as Fi = ff + A’‘e^ + (* = 1,2, ...,iV). 

Then we measure the variation of the potential energy from 
the minimum value, A$ = 'F — as a function of the am¬ 
plitudes, A^,Ao. If the mode k is harmonic. Ail? shows a 
parabolic shape. A"!* = [A^lo^Y/2, independent of A^, which 
is indicated by the black solid line. 


that occurs around t ~ 3. This time-scale coincides with 
the frequency of the lowest-lying mode, uj^ ~ 27r/t ~ 
2 , indicating that the increase in (Ar^(t)) is driven by 
large anharmonic contributions coming from the low-w^ 
modes. 



ji^k^k 



jik^k 



y^k^k 


FIG. 9. The potential energy landscape, A<I>(j 4*’, Ao), along 
(a) low = 2.19 and (b) high — 8.91 modes in crystal 1, 
and (c) low uj^ = 2.12 mode in crystal 2. See also the caption 
of Fig. [8l 


2. Potential energy landseape along eigen mode k 

In this section, we will interpret the above observation 
of strongly anharmonic vibrations (Fig. [5]) in terms of the 
local geometry of the energy landscape, i.e. variations in 
the potential energy compared to the parabolic, harmonic 
limit. In order to explore the energy landscape close to 
the minimum, we statically displace the particles from 
the energy minimum position along mode k, 

rl(A'=) = r1’+7l'=et (f = 1,2,..., iV). (17) 
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ri{A^) is therefore a function of the amplitude of 
mode k, and we vary A^ continuously from A^ = 0. 
In addition, to probe interactions between modes k and 
I ^ k, we can statically “excite” both modes and follow 
the displacements of particles that we can now express 
as 

( 3Af-3 ^ \ 

E ^ 

1=1,i^k ^ J 

The value of determines the extent of the excitations 
of modes I ^ k, and fi{A^,Ao) is then a function of both 
A’^ and Ao- Here it is worth to note that if modes I 
are harmonic in the T > 0 equilibrium state, their time- 
averaged amplitudes are = i/Tju;^ as in Eq. (IT4l) , thus 
the value of Ao provides a measure of the square root of 
the temperature, Ao = VT. The potential energy devi¬ 
ation, A<i> = $ — measured relative to the minimum 
value (presented in Table U), is then obtained as a 
function of the two amplitudes, A^, Ao; 

^ (19) 

(iJ) 

where rij{A^,Ao) = \ri{A'‘,Ao)-rj{A^,Ao)\, and 
is the summation over all the pairs of particles 

As in Eq. (EH), the potential energy landscape along 
the mode k is formulated in the harmonic approximation 

as fliiill; 

1 ^ ^ \0 
i—l j—\ ^ 

2 \df,dfj ^ 2 ’ 

i—l j — l ^ J / 

( 20 ) 

where Ui{A^,Ao) = fi{A^,Ao) — rf. Thus, if mode k is 
harmonic, (i) its energy landscape has a parabolic shape 
as a function of A^, and (ii) it is independent of Aq, 
i.e. it is independent of the other modes I ^ k. We 
plot Ad>(A^, Ao) for glasses in Fig.jSJ and for crystals in 
Fig. |9l where the harmonic behavior, A$ = (A^a;^)^/2, 
is indicated by the black solid line. 

Energy landscapes in glasses: Figure [5] shows the 
energy landscapes along, (a) low = 1.13 and (b) high 
uj^ = 8.79 modes in glass 1. The energy landscape of 
the low = 1.13 mode clearly deviates from a parabola 
shape. In addition, it is affected by the other vibrational 
modes I ^ k, depending on the value of Aq. We have con¬ 
firmed that the low = 1.13 mode is affected only by 
those low-lying modes with, ;< 5 i.e. mode couplings 
occur only between the lowest < 5 modes. In fact, 
the high = 8.79 mode maintains a parabolic shape re¬ 
gardless of the other mode excitations with Aq ^ 0. On 
the other hand, in glass 2, as seen in Fig. [SDc), even the 
low u!^ = 0.79 mode exhibits a harmonic energy land¬ 
scape, independent of Aq. These observations in glasses 


are consistent with the results of vibrational amplitudes 
reported in Fig. |nKa)-(c). The discontinuity in the inter¬ 
particle force in glass I (and also 3) deforms the energy 
landscape from a parabolic shape, and induces mode cou¬ 
plings. 

Energy landscapes in crystals: In Fig. [5J we plot 
the energy landscapes for crystals. In crystal 1, the en¬ 
ergy landscape along the low = 2.19 mode is rather 
flat and also sensitive to the other mode excitations I ^ k, 
resulting in strongly anharmonic behavior. The high 
= 8.91 mode shows a parabola shape for Aq = 0, but 
it is affected by the other mode excitations with Ao ^ 0. 
On the other hand, in crystal 2 (and also 3), all the 
modes become harmonic with a parabolic energy land¬ 
scape. This result is also consistent with the observation 
in Fig. inKd)-(f). 

One remarkable observation in crystal 1 is that the 
mode coupling due to cut-off nonlinearities are found 
even between low and high modes. This is not the 
case in glasses 1 and 3, where mode couplings occur only 
between the low modes. The interaction between the 
modes k and I emerges when those two modes share pairs 
of particles which experience the force discontinuity at 
r = Tc. In crystal 1, if two modes k,l have the same 
polarization and propagating direction, they share such 
pairs of particles, regardless of their frequencies, and 
On the contrary, in glasses 1 and 3, high modes 
are not extended, consisting of short spatial correlations 
(see Fig. [5]), and therefore such modes are less likely to 
share particles that experience the force discontinuity. 
Whereas, it is more likely for the low modes, which 
are extended in character, to have some overlap in their 
particle vibrations that experience the force discontinu¬ 
ity. 

In conclusion to this section, potential 2 (continuous 
interparticle force) exhibits no visible anharmonicities in 
the vibrational amplitudes as well as in the energy land¬ 
scapes along the modes. In contrast, there are clear 
signatures of anharmonic vibrational properties for po¬ 
tentials 1 and 3 (discontinuous interparticle forces). We 
emphasize again that for potential 2, although the in¬ 
terparticle force is continuous at r = Tc, cut-off non- 
linearities are induced by a second derivative which is 
non-analytic in r = r^. Hence, our results demonstrate 
that it is primarily the discontinuity in the interparti¬ 
cle force that enhances the cut-off nonlinearities, deforms 
the energy landscapes of many modes away from the ex¬ 
pected parabola shape, causing leaks in the energy from 
one mode to other modes, and as a result, induces de¬ 
tectable anharmonic effects. The number of particles 
pairs that experience the force discontinuity determines 
the strength of the cut-off nonlinearities for each mode 
k; more pairs of particles pass through the cut-off point 
r = rc, stronger anharmonic effects appear. 
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FIG. 10. The RDF g{r) for (a) glass 1 and (b) crystal 1 at 
zero, T = 0, and finite, T > 0, temperatures. The inset to 
(b) is a close-up around r = Vc = 2.5. The data at T — 0 are 
same as those presented in Fig. [2] 


C. Effects of cut-off nonlinearities on physical 
quantities 

In this Section we study the effects of cut-off non- 
linearities on several physical quantities, including the 
RDF, effective vDOS (Fourier transform of the velocity 
auto-correlation function), total potential energy, specific 
heat, and elastic constants. We will see that anharmonic 
effects manifest more evidently in quantities which in¬ 
clude more detailed information on vibrational states. 




FIG. 11. The first moments, uii and uji, of the vDOS, 
d{Ljj) and d{Lj), for (a) glasses and (b) crystals. We plot the 
temperature T > 0 dependence of wi, by the symbols, and 
uji, i.e. the values in the harmonic limit T = 0, by the lines. 

with that at T = 0. In crystal I, we observe the single 
peak around each lattice position at higher T > 10“^, 
while this peak splits into two peaks at lower temper¬ 
ature, T < 10“^ (see inset). These twin peaks are not 
related to the cut-off nonlinearities, but they simply orig¬ 
inate from the instability due to the force discontinuity at 
r = Tc, as discussed in Fig. [2Kb) and Sec. IIII A II There¬ 
fore, cut-off nonlinearities (see Figs. initoini) have no effects 
on g{r). We draw a similar conclusion for potentials 2 
and 3. 


1. Radial distribution function (static structure) 

In Fig. dQi we compare the g{r) a.t T > 0 with that 
at T = 0, for potential 1 (glass 1 and crystal 1). As the 
temperature increases, particles execute larger thermal 
vibrations, and the peaks in g(r) decrease in height and 
broaden. This temperature effect is clearly observed in 
the crystal case (Fig. fTUI b')'). There is also a decrease 
of the first peak in glass 1, although this change is quite 
small compared to the former case fFig. fTUI a')'). Figure [TUI 
demonstrates that as T tends to zero, g{r) smoothly con¬ 
verges to that at T = 0; g{r) at T = 10“^ coincides well 


2. Vibrational density of states 


Next we look at the effective vDOS, d{uj), obtained 
as the Fourier transform, of the velocity auto-correlation 
function d(t) at T > 0; 


d{t) 

dioj) 


N 




^(fi)(t)-u)(0)), 


i=l 



d(t) eyi'p{jujt)dt, 


( 21 ) 
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where Vi{t) is the velocity of particle i. For T —>• 0, d{uj) 
is expected to converge to the harmonic vDOS, d(a;), ob¬ 
tained from the dynamical matrix, Eq. (O [28| . This is 
tested in Figs. 151 and ITT] Figure [3] compares d{uj) (sym¬ 
bols) and d{uj) at T = 10“^ (lines), for the three investi¬ 
gated potentials. In addition, in Fig. [Til we compare the 
first moments, wi (lines) and uji (symbols), for d{uj) and 
J(w), respectively 

f'OO POO 

uJi= ujd{uj)duj, = / u}d{uj)duj. (22) 

do do 

In Fig. fTTT bl. we see small differences between d{ijj) and 
d{uj) for crystal 1; the value uii at finite T > 0 is slightly 
lower than wi at T = 0. This is caused by the anhar- 
monicities observed in Figs. [Bid) and|nia),(b). Although 
crystal 1 shows large cut-off nonlinearities for many eigen 
modes, they have only a small effect on the vDOS, shift¬ 
ing the mode frequencies to slightly lower values. Except 
for these small differences in crystal 1, Figs. 131 and [TT] 
demonstrate good agreement between d^uj) and d^uj) over 
the temperature regime studied. Therefore, we conclude 
that the vDOS and the actual values of the mode fre¬ 
quencies uj^ are insensitive to the cut-off nonlinearities. 


3. Potential energy and specific heat 




We next study the total potential energy deviation 
(A$) (measured from the minimum value $*^) and the 
specihc heat Cy, which are calculated as an ensemble 
average over the (NVT) MD trajectories as 

(A$)=/ ^ ^(r,,)\-$o, 

/ (23) 

{{E-{E)f) (SE^) 

'^V rp2 rp2 ’ 

where E = K [K = "^1/2 is the kinetic energy) 

is the total energy of the system, and the 5E = E — (E) 
are the associated fluctuations. 

In Figs.|T2]and[T3]we show the T-dependences of (A$) 
and Cv, respectively. Equipartition of energy directly 
provides the harmonic values for (A<I>) = (3iV — 3)T/2 ~ 
(3/2)iVT and Cy = 3A — 3 ~ 3A, which are indicated 
by horizontal solid lines in the figures. We clearly see 
that (A<h) and Cy coincide with these values in the in¬ 
vestigated T-range, for all the glasses and crystals, which 
is not surprising. Indeed, the potential energy (and the 
specific heat) is related to the eigen frequencies uj^ in the 
harmonic limit (Eq. (IT^ l. and the values of and g{ui) 
are not affected by cut-off anharmonicities as shown in 
Figs. 131 and nn 

It is, however, still surprising that crystal I retains 
strongly harmonic character in the values of (A$) and 
Cy, despite strong anharmonic effects which largely de¬ 
form the energy landscapes of the vibrational modes 


FIG. 12. The temperature, T, dependence of the scaled 
potential energy, (A<1?) /NT, for (a) glasses and (b) crystals. 
(A<1?) is the total potential energy measured from the mini¬ 
mum value (see Eq. (I23II 1. If thermal vibrations are har¬ 
monic, (A<1?} /NT = 3/2 - shown as the horizontal solid line. 


as shown in Figs. |9l(a) and (b). To further elucidate 
this point, we calculated the potential energy by using 
“parabolic” and “no-coupling” formulations; 


A$parabolic(t) = ^ ^-, 

k=l 

3N-3 

A$„o-coupling(t) = A$(A'=(t), Ao = 0), 

/c=l 


(24) 


where the value of A^(f) is the amplitude of mode k de¬ 
fined in Eq. ([T0|) . and A$(A*^,Ao) is the energy land¬ 
scape defined in Eq. ([T31). We recall that, as shown in 
Figs. l8l and [9] ISec. IIII B 2|) . cut-off nonlinearities (i) de¬ 
form the energy landscape from a parabolic shape and (ii) 
induce couplings between the modes. A$paraboiic disre¬ 
gards both effects, while A$no-coupiing includes (i) but 
not (ii). 

In Fig. [TTl we compare the time dependencies of 
those two quantities (A4>paraboia(t) and A<I>no-coupiing(t)) 
with that of the true value evaluated directly from the 
MD trajectory We see that A$paraboia(0 and 

A4>no-coupiing(t) Coincide well with the true A<l)(t) in glass 
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FIG. 13. The temperature, T, dependence of the specihc 
heat per particle, Cy jN (see Eq. (1231) 1. for (a) glasses and 
(b) crystals. If thermal vibrations are harmonic, Cy jN = 3 - 
shown as the horizontal solid line. 


1 (Fig. [Te a'll. In glass 1, even A<I’paraboia(i) works well to 
catch the instantaneous value of A<i>(t), because anhar- 
monic vibrations appear only in the low < 5 regime, 
and their effects are averaged out by the summation over 
all the modes (fc = 1 to 3A — 3). 

On the other hand, in crystal 1, the parabolic for¬ 
mulation A4>paraboia(t) Completely misses the true value 
A$(t) (Fig. [T^ blb This is obviously due to the large 
anharmonic amplitudes, A^{t). The no-coupling formu¬ 
lation A<I’no-coupiing(i), in contrast, behaves much bet¬ 
ter, since it takes into account the flattening of the 
energy landscapes due to large vibrational amplitudes, 
but it still deviates from the true value A<i>(t). Thus, 
both the flattening of the energy landscapes and the 
mode-couplings contribute to the total potential energy 
A$(t). As a final observation, as demonstrated by the 
MD trajectory value in Fig. im bl (red solid line), A$(t) 
fluctuates around the harmonic value (3/2)AT. Al¬ 
though the dynamics projected on the eigen modes de¬ 
fined at zero temperature does not follow the equipar- 
tition of energy, these low temperatures are enough 
to actually smooth out the energy landscape and re¬ 
sult in an effectively harmonic behavior for the en¬ 
ergy, in the spirit of the self-consistent phonon the- 
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FIG. 14. Time dependence of the potential energy per par¬ 
ticle, A<I>(t)/A, for (a) glass 1, (b) crystal 1, and (c) crystal 
2, at T = 10“®. Comparison between the true value from 
the MD trajectory (red solid line), parabola value A<I>paraboia 
(green dashed line), and no-coupling value A4>nc>.coupiing (blue 
dotted line) (see Eq. (1241) 1. Data are plotted every time in¬ 
terval At = 10. The horizontal black solid line indicates the 
harmonic value, A<F/A = {3/2)T = 1.5 x 10“®. 


ory [ 3 ^. Finally, as can be expected, crystal 2 shows 
good agreement between the three time histories of 

A$(t), A$parabola(t), Ad> 

no-cou pling(t) (Fig.EKc)). 
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FIG. 15. Dependence of non-affine moduli (Eq. (IA3ll 1 on temperature; (a) glass 1, (b) glass 2, (c) glass 3, (d) crystal 1, 

(e) crystal 2, and (f) crystal 3. We plot the average (global) value and the standard deviation 5M^ of the distribution 
function (see Eq. (IA2I) ~1. Values for the bulk (red circle), pure shear Gp (green square), and simple shear 

(blue triangle) moduli are presented. In glasses (a)-(c), the values for the two shear moduli, Gp and G^, coincide with each 
other, which is denoted by G^ (green square). The lines indicate values in the harmonic approximation, T = 0 (Eq. (IA4ll 1. 
Details of the calculations of elastic moduli are given in Appendix [Al 


4- Global and local elastic moduli 

As we have demonstrated above, effects of the cut¬ 
off nonlinearities are not particularly noticeable on the 


RDF, vDOS, total potential energy, and specific heat. In 
contrast, we show here that the elastic constants are rel¬ 
atively sensitive to cut-off nonlinearities. We computed 
the elastic modidi at finite T > 0, and at T = 0 in 
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the harmonic approximation, by using the methods de¬ 
scribed in Appendix In the present study, we focus 
on the non-affine components, and 

G^™ of the bulk, pure shear, and simple shear moduli, 
respectively, since anharmonic vibrations are directly re¬ 
flected in these non-affine components. Here m indicates 
coarse-grained cubic domains of size much smaller than 
the simulation box size, as detailed in Appendix The 
affine components, , G^ , G^™, are deter¬ 

mined rather by static structural properties, which are 
insensitive to anharmonicities as seen in g{r) of Fig. [TUI 
We have indeed confirmed that the affine components are 
not influenced by cut-off nonlinearities. 

In Figure [TU| we present the average (global) value 
and the fluctuations (standard deviation) 5M^ of 
the probability distributions (see Eq. (IA2l) b 

as functions of T, for each considered glass and crystal. 
5M^ is a measure of the elastic heterogeneity. In the 
glasses (Fig.|TUDa)-(c)), both values of and 6M^ are 
insensitive to T variations over the temperature regime 
studied. Glass 1 clearly shows a disagreement between 
the finite T > 0 value and the harmonic value (T = 0), 
whereas good agreement is observed in glass 2. The dis¬ 
agreement in glass 1 is caused by the cut-off nonlinear¬ 
ities, which are enhanced by the discontinuous force at 
r = rc- We note that the difference between the T > 0 
and T = 0 values shows up more clearly in the shear mod¬ 
ulus G^,6G^ than the bulk modulus K^,SK^, since 
the shear modulus expresses a larger non-affine compo¬ 
nent In addition, 6M^ shows larger discrepan¬ 

cies than the global , indicating that local quantities 
are more sensitive to anharmonicities than global values. 
The cut-off nonlinearities cause the T > 0 value of SM^ 
to increase away from the T = 0 value, i.e. the moduli 
distributions become more heterogeneous. Indeed, in dis¬ 
ordered amorphous structures, particles feel the force dis¬ 
continuity at r = Tc randomly in space. Cut-off nonlin¬ 
earities therefore induce spatially heterogeneous effects. 

Glass 3 also shows disagreement between the T = 0 
and T > 0 values of ,SM^, albeit a smaller difference 
than that for glass 1. In glass 3, cut-off nonlinearities are 
smaller, and have a smaller effect on the elastic moduli. 
In both glasses 1 and 3, the disagreement persists down 
to the lowest temperature, T = 10“^, which is three or¬ 
ders of magnitude lower than Tg ~ 0.4. This observation 
therefore poses a caveat on the value T = 10“^, which 
is sometimes considered as a temperature where the har¬ 
monic limit is recovered. As we have shown above, this 
peculiar behavior is due to the fact that cut-off nonlinear¬ 
ities persist to the very-low-T regime, contrary to usual 
expansion nonlinearities which vanish at small tempera¬ 
tures. 

Turning to the crystals iFig. fTSl dVff')'). we realize that 
the cut-off nonlinearities produce completely unexpected 
and anomalously large values of the elastic moduli in 
crystal 1. At low-T, the mechanical response of crys¬ 
tals is described only in terms of the affine modulus |35l | . 
Additionally, the ordered lattice structure of a uniform 


crystal should lead to a spatially distribution in the elas¬ 
tic moduli that is homogeneous. As a result, we expect 
that not only should the (total) non-affine modulus van¬ 
ish, = 0, but so should the corresponding fluctu¬ 
ations, SM^ = 0, i.e. spatially, all the local values of 
non-affine moduli, = 0. In fact, this is indeed 

true for crystals 2 and 3, as demonstrated in Fig. me) 
and (f). Note that the increase of at T = 10 ^ is 
caused by expansion nonlinearities due to large thermal 
vibrations. However, crystal 1 clearly shows large devi¬ 
ations in > 0 and 6M^ > 0 compared to the zero 
expected value. For the case of the shear moduli fluctu¬ 
ations, we note that although the values are small, SGp 
and SGf assume finite values even at T = 0. This comes 
from the fact that particles are slightly displaced from 
the exact lattice positions due to the force discontinuity 
at r = Tc = 2.5, i.e. crystal 1 is slightly inhomogeneous, 
as demonstrated by the g{r) of Fig.jUKb). 

We underline here an interesting observation in the 
case of crystal 1; the bulk and pure shear G^ mod¬ 
uli are strongly affected by the cut-off nonlinearities, 
whereas the simple shear modulus Gf is much less sensi¬ 
tive. This result indicates that the simple shear phonon 
modes, i.e. transverse phonons in the direction [1,0,0], 
tend to avoid the force discontinuity at r = Tc and are 
less affected by the cut-off nonlinearities, compared to 
phonons propagating in other directions. 

To sum up this section, anharmonicities due to cut-off 
nonlinearities cause relatively large effects on the elastic 
moduli. In the glassy states, systematic deviations up to 
10 -20% persist in the moduli and their fluctuations over 
the entire temperature range. While, for the crystals, 
anharmonic effects can cause significant deviations even 
at the lowest temperatures. None of these anharmonic 
effects are evident in the RDF, vDOS, potential energy, 
and specific heat. These facts point to the conclusion 
that the strength of the anharmonic effects on a physical 
quantity seems to be determined by the total amount of 
detailed information regarding vibrational excitations it 
includes. Indeed, the vDOS, potential energy, specific 
heat, only involve functions of the eigen frequencies 
and are only mildly sensitive to the cut-off nonlinearities. 
In contrast, as it is evident in the harmonic equation, 
Eq. dm, the non-affine moduli are described in terms 
of both the eigen frequencies and eigen vectors e^, 
i.e. they include a complete information on the structure 
of the vibrational modes. Also note that eigen vectors 
are more strongly affected by the nonlinearities than the 
eigen frequencies Q- 


IV. SUMMARY 

In the present paper, we have studied the effects of 
cut-off nonlinearities induced by the interaction potential 
cut-off on low-temperature thermal vibrations in glasses 
and crystals. It is common practice in computer simu¬ 
lations to truncate the (long-range) interaction potential 
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at some cut-off distance r = Vc- Here, we have focused 
on three specific parameterizations of the traditional LJ 
potential, as listed in Table U and Fig. m (i) potential 
1: shifted potential </>sp(r) with = 2.5, discontinu¬ 
ous 1st and 2nd derivatives, (ii) potential 2: shifted-force 
potential (j)SF{r) with rc = 2.5, discontinuous 2nd deriva¬ 
tive, and (iii) potential 3: shifted potential ipspir) with 
Tc = 3.0, discontinuous 1st and 2nd derivatives. While 
these three truncated potentials do not show any notice¬ 
able differences in the static structure (Fig. [5]) and vi¬ 
brational states (Figs. [3] and 2]) in the (harmonic) limit 
T = 0, differences become apparent in thermal vibra¬ 
tions at finite temperatures T > 0 (Figs. [6] to |9]), which 
originate from the cut-off nonlinearities. The truncation 
causes the potential to be non-analytic at the cut-off, and 
the harmonic equations of motions, Eq. o, do not ac¬ 
count for this singularity. Thus, when some pairs of par¬ 
ticles pass through the harmonic description breaks 
down, leading to anharmonicities in physical observables. 
The cut-off nonlinearities are distinct from the usual ex¬ 
pansion nonlinearities, which become more prominent at 
higher temperatures, that come from neglecting higher 
order terms in the Taylor expansion of the potential. 

We have demonstrated that it is primarily the discon¬ 
tinuity in the interparticle force (the first derivative of 
the potential) at r = rc that enhances the cut-off non- 
linearities. The force discontinuity deforms the potential 
energy landscapes of many eigen modes away from the 
harmonic parabolic shape, as shown in Figs. [5] and [51 In 
addition, when normal modes share the same pairs of par¬ 
ticles experiencing the force discontinuity, they exchange 
or leak energy, leading to further deformation of their en¬ 
ergy landscapes due to the modes coupling. These effects 
can be quantified in terms of the vibrational amplitudes 
shown in Fig. [51 For crystal 1, where the force discon¬ 
tinuity is located at = 2.5, which also corresponds to 
the nearest-neighbor spacing, the implication of this are 
particularly magnified precisely because the distribution 
of particle separations is a delta function at Tc itself. As 
a result, most, if not all, particle pairs participating in 
normal mode vibrations experience the force discontinu¬ 
ity at the same time, leading to the observed enhanced 
anharmonic effects. On the contrary, even though cut-off 
nonlinearities also exist in glass 2 and crystal 2, where the 
force is continuous at r^, anharmonic effects were largely 
suppressed and barely detectable in the vibrational am¬ 
plitudes and energy landscape features. 

We also studied the effects of cut-off nonlineari¬ 
ties on several physical quantities, including the RDF, 
vDOS, potential energy, specific heat, and elastic mod¬ 
uli (Figs. [10] to [151). We showed that the RDF, vDOS, 
potential energy, and specific heat are resilient to cut-off 
effects, whereas the elastic moduli are rather sensitive to 
the anharmonic nature of the normal modes. Indeed, the 
iron-affine components of the elastic moduli, computed 
within the harmonic approximation (T = 0), are deter¬ 
mined by both the eigen frequencies and crucially 
the eigen vectors e* (see Eq. (El), i.e. they include 


1.6 



FIG. 16. Amplitude versus eigen frequency for 
the shifted potential (j>sp{r) with different cut-off distances, 
Tc = 2.5, 3.0, 4.0, and 5.0. The configuration is glass, and 
the temperature is T = 10“®. Note that rc = 2.5 and 3.0 cor¬ 
respond to glasses 1 and 3, respectively. See also the caption 
of Fig. El 

complete information on the structure of the vibrational 
excitations. Therefore, it is understandable that anhar¬ 
monic effects have a relatively large impact on the elastic 
moduli. Particularly, the large anharmonic vibrations in 
crystal 1 results in unphysical trends in the elastic mod¬ 
uli. In contrast, the vDOS, potential energy, and specific 
heat, which include only the eigen frequency information 
in the harmonic limit, are insensitive to anharmonicities. 

In contrast to expansion nonlinearities, cut-off non- 
linearities persist in the very-low-temperature regime, 
T = 10“^ <C Tg, Tm- This is especially the case for 
crystal 1, where we have found that cut-off nonlinearities 
actually become more prominent with decreasing temper¬ 
ature. As a result, the elastic moduli of glasses 1,3 and 
crystal 1 do not coincide with the harmonic values even 
at T = 10“^, where we normally expect the harmonic 
description to be valid. Even though thermal displace¬ 
ments are very small, Ar ^ \/T, pairs of particles with 
r = rc feel the force discontinuity, causing the anhar¬ 
monicities. This counter-intuitive behavior of the cut-off 
nonlinearities is obviously at odds with our understand¬ 
ing of expansion nonlinearities, which vanish in the low 
temperature limit. 

V. CONCLUSIONS 

We conclude with a few observations. The factor ul¬ 
timately determining the strength of the cut-off non- 
linearities is the number of interacting particles that ex¬ 
perience the cut-off discontinuity or, more specifically, 
the number of inter-particles distances close to the cut¬ 
off distance. This in turn depends on the physics of the 
system under study. For the case of the (traditional) 
12-6 LJ potential explored here, we found strong dis¬ 
crepancies for crystal 1 because of the commensurability 
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between the cut-off distance and the choice of the lat¬ 
tice environment (that was tuned according to the den¬ 
sity) . The anharmonicities arise because particles fall out 
of contact with one another during thermal vibrations. 
Thus, from an interaction point of view the potential ap¬ 
pears to behave one-sided for particle pairs separated by 
the cut-off distance. This situation is similar to models 
of granular materials and jammed particulate packings, 
where particles interact thr ough one-sided, finite-range 
(contact/cut-off), potentials [2^- 

Also, for studies at low temperature, in or close to the 
harmonic regime, one has to be extremely careful with 
the implementation of the potential truncation. A cut-off 
distance Tc = 2.5cr is often employed in computer simu¬ 
lations involving LJ-type potentials flol - [l^ . However, as 
we have seen in the present work, this value may not be 
suitable to study low-temperature vibrational properties 
for particular systems, as it can result in inconsistencies 
between the harmonic value T = 0 and the low temper¬ 
ature values. Though we remark here that in the case 
of the soft-sphere potential, (jiscif) = which is 

popular repulsive potential employed in many numerical 
works [l3, 11^, the cut-off discontinuity is significantly 
smaller than for the 12-6 LJ potential. In this case, the 
value of (fisci'r) at = 2.5 is 1.7 x 10“^% of e, which 
is much smaller than that of the LJ potential, 1.6% of 
e. In our previous works [13,1131, we employed = 2.5 
in (f’sci'i’)-: which was large enough to study both elas¬ 
tic constants and sound and thermal transport, without 
any visible cut-off nonlinearities. Here, we have observed 
that physical quantities that depend only on the struc¬ 
tural arrangement of the particles and/or just the eigen 
frequencies, such as the RDF, vDOS, potential energy, 
specific heat, are not appreciably sensitive to the cut-off 
effects. However, cut-off nonlinearities may be apparent 
in others quantities that depend explicitly on the eigen 
vectors, or particle polarization vectors of the individ¬ 
ual particles, such as mechanical properties, for example 
the elastic moduli. Thus, the apparent validity of the 
harmonic r^ime can depend on which observable is con¬ 
sidered [a III. Low frequency modes are more influenced 
by the cut-off nonlinearities, so that thermal transport 
properties are also likely to be affected in a signihcant 
manner. 

Ideally, we would like to simulate a LJ system with¬ 
out any cut-off of the potential, or one with a long cut¬ 
off distance. Indeed, the cut-off nonlinearities can be 
suppressed by employing longer cut-offs, thereby reduc¬ 
ing the discontinuity at r = Cc. In Fig. [12] we plot 
the vibrational amplitudes at the indicated values 
of Tc = 2.5, 3.0, 4.0, and 5.0, for the glass configura¬ 
tion at temperature T = 10“^. From these data, we see 
that the values of Vc = 4.0 or 5.0 can be sufficient to 
suppress the anharmonic effects in the present LJ sys¬ 
tem. However, we have shown that it is the discontinuity 
in the first derivative of the potential that has a domi¬ 
nant influence. In addition, we have also confirmed that 
the vibrational modes are practically identical in the case 



T 

K 

Lrp ^p 

/^A /-iN 
\jr s Kjr g Kjr g 

Glass2 

0 

60.7 61.2 0.5 

14.0 35.7 21.7 

14.0 35.7 21.7 

10"® 

60.7 61.2 0.5 

13.9 35.7 21.8 

13.9 35.7 21.8 

Crystal2 

0 

43.4 43.4 0.0 

14.4 14.4 0.0 

38.1 38.1 0.0 

10"® 

43.7 43.7 0.0 

14.7 14.7 0.0 

38.2 38.2 0.0 


TABLE II. The average (global) elastic moduli, K, Gp, and 
Gs, of glass 2 and crystal 2. The temperature is T = 0 and 
10“®. We also present values of the affine and non-affine mod¬ 
uli. Note that for glass 2, Gp = Gs. 


of long cut-offs (re = 4.0, 5.0) and shifted-force potential, 
i.e. force-shifting procedure practically does not alter the 
vibrational modes. Therefore, it is recommended that the 
interparticle force (the first-derivative of the potential) be 
made continuous at r = Tc, as proposed in Refs. [ofl^. 
to avoid unwanted anharmonic effects without increasing 
the computational costs by using a large cut off distance. 

On the other side, for those potentials that are in¬ 
trinsically truncated (e.g. soft elastic particles), sys¬ 
tems are likely to exhibit nonlinearities im that can 
strongly influence mechanical properties in athermal or 
low-temperature states. Thence, it is particularly im¬ 
portant for such systems, that efforts, which extrapo¬ 
late transport and mechanical properties into the ther¬ 
mal regime purely from structural data or the static har¬ 
monic formulation, should be treated with caution. For 
instance, predicting material properties such as modes of 
failure could potentially lead to unexpected catastrophic 
effects when such results be employed in the design of 
critical devices. 


ACKNOWLEDGMENTS 

We thank H. R. Schober, M. Sperl, Th. Voigtmann, 
A. Ikeda, and C. S. O’Hern for helpful discussions and 
useful comments. 


Appendix A: Measuring finite temperature T > 0 
and zero temperature T — 0 elastic moduli 

We computed the elastic moduli at finite temperatures 
T > 0 as well as at zero temperature T = 0 (the harmonic 
limit), the results of which are presented in Sec. IHI (J4l 
In order to measure the elastic moduli, we did not apply 
any explicit deformation, but rather we employed the for¬ 
mulation developed from the linear response theory. For 
calculation of the finite T > 0 moduli, we used the equi¬ 
librium fluctuation formulae [ssl - l^ . which we referred 
to as “fully local approach” in our previous study [l8l| . 
For calculation of the zero T = 0 moduli, we extended 
the formulation for the global moduli, which has been 
established and used in previous works [43 - l43| . to the 
local moduli. In our recent work [4^, we have employed 
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this extended formulation to study local elastic moduli 
distributions in T = 0 athermal jammed solids. 

The present 3-dimensional system was subdivided into 
a grid 40 x 40 x 40 cubic domains of linear size w = 
3.16 (coarse graining length). For each cubic domain 
m {m = 1, 2,..., 64000), we measured the local modulus 
tensor, (a, /?, 7, <5 = x, y, z), which is defined as the 

derivative of the local stress tensor with respect to the 
linear strain tensor [l^ . The value of is composed 

of four terms; the Born term the kinetic term 

C^^g, the stress correction term C^^g, and the non- 
affine term C^^g] 


glasses, so we represent G™ for the shear moduli in 
glasses, i.e. P(G'") = P(G™) = P(G™). For ref¬ 
erence, Table El presents values of the global moduli, 
M = K,Gp,Gs, for glass 2 and crystal 2 at zero tem¬ 
perature, T = 0, and a finite temperature, T = 10“^. 

In the present study, our major focus is on non-alRne 
contributions G^^g, since anharmonic effects are directly 
reflected in G^^g, while affine components G^^g are 
mainly determined by static structural properties and are 
therefore rather insensitive to anharmonicities. At finite 
T > 0, G^^g is measured in terms of a local and global 
stress correlation function [3§ - l4lj| : 


C m _ /^Bm I /^Km , /^Cm ^Nm 

apyS ^O'13 'yS ' ^a.0'y5'> / A 1 ^ 

C Am /^Nm '' ' 

q ;/37(5 ^a.^'^5' 

The summation of the first three terms, i.e. G^^g = 

is the “afRne” modulus [13, 
[37I . [4]] |. Here we note that the Born term G^^g is de¬ 
fined as the second derivative of the energy density with 
respect to the Green-Lagrangian strain tensor [3^ [43 . 
Therefore, if we define the modulus by using the linear 
strain tensor, the stress correction term G^^g is nec¬ 
essary as long as the stress tensor has finite values in 
its components [i^. The affine modulus describes the 
elastic response, when particles follow the applied affine 
strain field and are displaced affinely at all scales. On 
the other hand, the forth term, G^^g, is referred to as 
the “non-affine” modulus, which contributes negatively 
to the overall modulus, and comes from “additional” par¬ 
ticle displacements at the microscopic scale that deviate 
from the applied affine field [^. Finally, from the 
modulus tensor G'^p^g, we calculated the bulk modulus 
AT"*, pure shear modulus G™, and simple shear modulus 
G^ in the same way as in Refs. [H, [13,1131 • 

After measuring the local elastic moduli in each little 
cube m (to = 1,2,..., 64000), we collected data to obtain 
probability distribution functions P(M™) for the mod¬ 
uli, M"* = A:™,G;7,G™. From P(M™), we calculated 
the average (global) value M and the standard deviation 
(fluctuations) 5M] 


M = 




SM = 





(A2) 


Two shear moduli distributions, P{G^) and P{G^), 
coincide with each other in the isotropic systems like 


CXs = ^[«P<^^s)-Kp){<J,s)], 
= 7 {S<ySa,s) , 


(A3) 


where Uap and are the global and local stress tensors, 
respectively (see Ref. [l^ for a detailed formulation), and 
5(Ja.p = (Xa.p - {(Tap) and - {(Xap^ are the 

corresponding fluctuations. () denotes the NVT ensem¬ 
ble average, which coincides with the average over the 
MD trajectory when the dynamics is ergodic 0- Since 
the function of G^^g contains three- and four-point cor¬ 
relations [13; El; we therefore need a long simulation 
trajectory for good numerical convergence of the ensem¬ 
ble average [l8| . To that end, we performed production 
runs for 2 x 10^ to 4 x 10^ simulation steps at each tem¬ 
perature r, and the averaging was performed over 10^ 
to 2 X configurations separated by 2 x 10^ steps for 
glasses, and over 10® to 2 x 10® configurations separated 
by 2 X 10^ steps for crystals. 

In the T = 0, harmonic limit, the non-affine term 
G^^g is formulated in terms of the eigen frequency 
and the eigen vector e* [i^ - fil] : 


/^Nm 


3N-3 

E 






d<p\ 

On ) 




OiV —o J ^ 

= E 

fc=l ^ 

(A4) 

where da^^p = ^Zi ^ ' daapidfi and da)^^ = Y{Zi ^ ’ 
da'^pldn are the fluctuations of the global and local 
stress tensors, induced by the eigen mode k. We note 
that Eq. (jA4p is obtained by taking the limit of T —0 
in the finite T > 0 formulation, Eq. (IA3I) [ 13 ]. 
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